function dist = catchSSBF_dist(est0,H0,SSB0,normH,normSSB,F_at0,normF,...
    R,T,w_a,m_a,mature)
% input: 
% est0 = initial guess of parameter values
% H0 = observed catch data (standerizded time series, kg)
% SSB0 = observed SSB (standerizded time series, kg)
% R = observed recruitment data (time series, number of fish)
% T = number of age classes
% w_a = weight per fish (kg)
% mature = dummy variable indicating age at maturity
% m_a = natural mortality

ttmax = length(H0); % number of years
A = size(F_at0,2);
age = 1:T;
m_a = m_a(1:T);
t1_range = (SSB0 >= 0); % years >= 1968
t2_range = (R <= 0);

% parameters to be estimated ==========================
R_scale = est0(1);
lambda = est0(2);
ft = est0(3:(ttmax+2));
Theta = est0((ttmax+3):end);
K = length(Theta);

R0 = max(R)*R_scale; % Exog recruitment before 1968
N0 = R0*exp(-lambda*age);

% compute age-dependent fishing mortality F_hat(t,a)
XX_age = [ones(T,1) age' (age.^2)' (age.^3)' (age.^4)' (age.^5)' (age.^6)' ...
    (age.^7)' (age.^8)' (age.^9)' (age.^10)'];
X_age = XX_age(:,1:K);
h_age = X_age*Theta; % size = [T 1]
h_age(A+1:end) = h_age(A); % set constant from age 10

% compute fishing mortality *******************************
F_hat = ft*(h_age');

z = F_hat + ones(ttmax,1)*m_a'; % total mortality (t,a)

% population matrix simulation ---------%
N = zeros(ttmax,T+1); % row= year, col=age
N(t2_range,1) = R0; % assume Recruitment to be constant and exog
N(t1_range,1) = R(t1_range);
N(1,2:end) = N0;
for t = 2:ttmax
    for j = 2:(T+1)
        N(t,j) = N(t-1,j-1)*exp(-z(t-1,j-1)); 
    end % for j
end % for t

% predicted catch & SSB ----------------------%
H_age = zeros(ttmax, T); % age-structured catch H(t,a)
H_age = (F_hat./z).*(1-exp(-z)).*(ones(ttmax,1)*w_a).*N(:,2:end);
H_hat = H_age*ones(T,1);


SSB_hat = zeros(ttmax,1); % predicted spawning stock biomass
SSB_hat = N(:,2:end)*(mature.*(w_a'));

v1 = H0 -(H_hat./normH); 
v2 = SSB0(t1_range)-(SSB_hat(t1_range)./normSSB);
v3 = F_at0(t1_range,1:A) - (F_hat(t1_range,1:A)./normF);
v3 = v3(:);

dist = norm(v1)+ norm(v2) + norm(v3); 
end